pro stochastic_sed_mktable

lsun=3.85d33
dir='$SLUG_DIR/data/020/final/'
tracks=mrdfits(dir+'tracks_020.fits', 1)
lej=mrdfits('lej_p00_final.fits', 1)
lej=lej[where(lej.teff GT 0)]

log_g_list=unique(lej.logg, /fin)
log_teff_list=alog10(unique(lej.teff, /fin))
nteff=n_elements(log_teff_list)
ng=n_elements(log_g_list)
nflux=n_elements(lej[0].flux)

table=fltarr(nteff, ng, nflux)
table_bol=fltarr(nteff, ng)
for i=0, nteff-1 do begin ;populate table
  lejwh=lej[where(alog10(lej.teff) EQ log_teff_list[i], ct1)]
  if ct1 NE 1 then begin
   for j=0, nflux-1 do table[i, 0:ng-1, j]=interpol(lejwh.flux[j], lejwh.logg, log_g_list)
   table_bol[i, 0:ng-1]=interpol(lejwh.bolometric_lum, lejwh.logg, log_g_list)
  endif else begin
     table[i, 0:ng-1, 0:nflux-1]=lejwh.flux##replicate(1, ng)
     table_bol[i, 0:ng-1]=lejwh.bolometric_lum
  endelse
endfor
STOP
big_table=fltarr(n_elements(tracks), n_elements(tracks[0].age), n_elements(lej[0].flux))


for i=0, n_elements(tracks)-1 do begin
  for j=0, n_elements(tracks[0].log_g)-1 do begin
	xind=interpol(findgen(nteff),log_teff_list, tracks[i].log_teff[j])
	yind=interpol(findgen(ng),log_g_list, tracks[i].log_g[j])
	lum00=interpolate(table_bol, xind, yind)
	sed=fltarr(nflux)
	for k=0, nflux-1 do sed[k]=interpolate(table[0:nteff-1,0:ng-1,k], xind, yind)
	scl=lsun*(10.^(tracks[i].log_lum[j]))/lum00 ;scale by the luminosity
	big_table[i, j, *]=sed*scl
   endfor 
endfor


time_list=unique([tracks.age])
ntime=n_elements(time_list)
mass_list=reverse(unique(tracks.mass_star));because tracks are opposite direction
nmass=n_elements(mass_list)
table_final=dblarr(ntime, nmass, nflux)
ind_end_ms=50

for i=0., ntime - 1 do begin
  t=time_list[i]
  for j=0, nmass-1 do begin
    if max(tracks[j].age[ind_end_ms]) GT t then begin
        for k=0, nflux-1 do table_final[i, j, k]=interpol(big_table[j,*,k],tracks[j].age,t) 
    endif  else table[j, *, *]=0
  endfor
  print, i/(ntime-1.)
endfor

STOP
mwrfits, table_final, 'sed_table_final.fits', /create


end
